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Abstract. The weighted density approximation, its implementation and its applica- 
tion to ferroelectric materials is discussed. Calculations are presented for several per- 
ovskite oxides and related materials. In general the weighted density approximation is 
found to be superior to either the local density or generalized gradient approximation 
for the ground state. Electronic structures are little changed. The linear response of 
the weighted density approximation is calculated for the homogeneous electron gas, and 
found to be improved relative to the local density result, but not in full agreement with 
existing Monte Carlo data. It is shown that the agreement can be further improved by 
a simple modification. Calculations of the ferroelectric soft mode in KNbC>3 suggest 
that the low temperature distortion is approximately 20% smaller than indicated by 
existing experiments. 



I INTRODUCTION 

Piezoelectric, ferroelectric and related perovskites are both technologically im- 
portant and physically interesting systems and, as such, have been the subject of 
considerable recent interest. Several groups have performed first principles calcu- 
lations based on density functional theory. These calculations, primarily based on 
the local density approximation (LDA) have elucidated important aspects of the 
underlying physics of these materials as well as providing quantitative information 
about phonons, polarizations, crystal structure, elastic and other properties. [1] 

Two of the many themes that have emerged are (1) Ferroelectric and piezoelec- 
tric properties are due to a delicate balance between large competing interactions, 
[2,3] so very high accuracy is required for predictive calculations and (2) The lat- 
tice instabilities behind these phenomena are extraordinarily sensitive to volume. 
In fact, in important materials, like KNbO"3 [4-10] and BaTiO"3 [2,3,11] ferroelec- 
tricity, while predicted at the experimental volume, is absent, or at least strongly 
suppressed, at the predicted LDA volume. 

Generally excellent agreement with experimental knowledge has been obtained 
when calculations are performed at the experimental volume. This includes even 
transition temperatures, polarizations and transition temperatures [1]. In anti- 
ferroelectric PbZrC>3, LDA calculations [12] at the experimental volume yielded a 
crystal structure that differed from existing experimental measurements. However, 



the LDA structure has now been confirmed by two independent neutron inves- 
tigations, [13,14] providing yet another demonstration of the utility of the first 
principles approach. On the other hand, even at the experimental volume, the 
LDA consistently overestimates the dielectric constant, 6^. [15-17] Exploration of 
novel systems that are of necessity less well studied experimentally than current 
materials will require an ability to calculate properties independent of experimen- 
tal measurements of the volume. This may also be important for determination of 
piezoelectric properties in materials like PbTi03 where the LDA gives an incorrect 
prediction of the strain associated with the deviation of the c/a ratio from unity 
even with the volume constrained. [18] This raises the question of whether there is a 
practical way of correcting the LDA errors, in order to construct a more predictive 
first principles approach for these materials. 

The LDA errors can sometimes be trivially fixed. For example, if the experi- 
mental volume is known, the calculations may be done constrained to this volume. 
In some cases, e.g. KNb0 3 , simple modification of the LDA functional through 
use of the Wigner form rather than an electron gas-like form, [7,8] or the use of 
generalized gradient approximations (GGA) to density functional theory are suf- 
ficient to correct the volume (although the dielectric constant remains too high). 
However, in general these approaches do not solve the LDA problem. For example, 
in BaTiO"3 the GGA overcorrects the LDA volume, so that the error is as large but 
in the opposite direction. [9] The implication is that a more sophisticated, presum- 
ably non-local, density functional is needed to fully circumvent the problems noted 
above with the LDA. 

The first efforts at developing practical non-local functionals date from the 1970's 
when the average density approximation [19] (ADA) and weighted density approx- 
imation [20-22] (WDA) were proposed. However, over most of the intervening 
period the field has been relatively dormant, in part because of the success of 
the simpler LDA and GGA schemes and in part because it was widely thought 
that such schemes could not be implemented in a computationally tractable fash- 
ion. However, at least for the WDA, computationally efficient algorithms are now 
known [10,23-25] and benchmark calculations have been reported. The ADA has 
attracted less interest. 

In the cases that have been studied ground state properties of solids are gener- 
ally improved over the LDA. [26] These include tests for several simple elements 
and compounds and KNb03. Significantly, this latter test showed that the WDA 
predicts an equilibrium volume for KNbOs in almost perfect agreement with exper- 
iment. [10] Unlike the GGA, which also greatly improves the equilibrium volume 
of KNbOs, the WDA is a truly non-local density functional, in that the exchange 
correlation potential at a point r incorporates information about the charge density 
n(r') over a finite region of space. 

Here, following a brief overview of the method, we report investigations of a 
range of perovskite and related oxides within the WDA. These show much im- 
proved volumes relative to either the LDA or GGA. Phonon frequencies and the 
ferroelectric soft mode in KNb0 3 are calculated within the WDA and compared 



with LDA results and experiment. An analytical expression for the linear response 
of the uniform electron gas in the WDA is derived. 

II THE WEIGHTED DENSITY APPROXIMATION 

The WDA and the ADA are both based on the general expression for the ex- 
change correlation energy of a general electron gas E xc , in density functional theoryi 
(DFT). 

E xc = y / ^^G(r, r'){n(r)}rfrdr', (1) 

where the function G(r, r') is also a functional of the total electronic density n(r). A 
rigorous expression for G can be derived [27] in terms of coupling constant averaged 
pair correlation function: 

G(r,r') = [\g( r , r ';\){n(r)}-l]d\. (2) 
Jo 

Were the actual pair correlation function used instead of the coupling constant aver- 
aged one, this expression would give the interaction energy of each electron with its 
exchange correlation hole. The difference due to the averaging reflects the fact that 
kinetic energy of interacting electrons differs from that of a non-interacting system 
with the same density. It is this difference, also known as exchange-correlation 
kinetic energy, which is responsible for the high q, short wavelength behavior of the 
response, and is also implicated in usual underestimation of semiconducting gaps 
in Kohn-Sham DFT. 

For the uniform gas this function, Gq(\t — r'\,n), is known with high accuracy 
[28] , but for an arbitrary inhomogeneous system, like occurs in a real density func- 
tional calculation for a solid, G is not known and there is therefore no practical 
way to use this formula without making some approximation. 

The LDA instead of Eq. (1) uses (e 2 /2) / drdr'n 2 (r)G [\r - r'|, n(r)]/|r - r'|, 
so that E xc becomes E^ A = J n(r)e xc [n(r)]dr, e xc being the density of exchange- 
correlation energy of the uniform gas. The LDA is incorrect in both most important 
limits: the fully localized, i.e., a one electron system, and the fully delocalized limit, 
i.e., homogeneous electron gas. In the former case the LDA gives a spurious self- 
interaction with energy (e 2 /2) / drdr'n(r)n(r')/\r — r'| + / n(r)e xc [n(r)]dr, which is 
widely thought [29] to be a key problem with the LDA. In the homogeneous limit, 
the LDA gives the correct exchange-correlation energy, but the changes of this 
energy upon small perturbations are not properly described; the second variation 
of E xc with density, i.e., the exchange-correlation part of the dielectric response, 
K xc (r — r') = S 2 E xc /Sn(r)Sn(r'), is a delta function, which is incorrect. The Fourier 
transform of K xc (r) in LDA is independent of the wave vector. Since LDA is exact 
for the uniform gas, K^ A corresponds to the correct K xc at q — 0. GGAs also give 
correct behavior at q — 0, but become even worse than the LDA at high g's. 



The two nonlocal expressions for E xc , WDA and ADA, are aimed at correcting 
one or the other of these two limits. The former uses the general expression (1), 
but instead of the actual function G uses a model function, defined so that the 
one electron limit is honored. This begins by choosing a generic expression for G, 
which depends on one parameter n, to be defined later. In the original papers it 
was suggested that G(r, r', n) = G (r, r', n) = fo[g(r, r'; A, n) — l]d\, where g is the 
averaged pair correlation function of the homogeneous electron gas. Later it was 
realized [30] that other choices of G may be better than Go, and that there is no 
physical reason to prefer Go over many other choices. In the WDA n is a function 
of r, but differs from n(r), and is chosen so that / G[r, r', n(r)]n(r')dr' = —1. This 
assures that for a one electron system E xc cancels the self-interaction exactly. 

Within the WDA, non-local information about the charge density is incorporated 
into E xc both through the construction of n via the sum rule, and through the non- 
local Coulomb integral. 

We note that G need not neccesarily be the actual pair correlation function of 
the system: although Eq.(l) has the same functional form as the WDA energy, 
the fact that G WDA is a function of an averaged density n, not a functional of 
the true density n(r), means that it is possible that the best approximations for 
qWda cou y kg different from the physical function G defined in Eq. (2) even for 
the uniform electron gas. However, the few reported calculations suggest that this 
function, while perhaps not optimal, does yield results much superior to the LDA. 
Besides, use of the homogeneous electron gas averaged pair correlation function has 
an appealing simplicity and calculations of ground state properties, below, use this 
ansatz. The results confirm that at least for the structural properties considered, 
the WDA with this choice of G yields results that are uniformly superior to the 
LDA. 

Ill GROUND STATE PROPERTIES IN THE WDA 

As mentioned, one of the basic problems with application of the LDA to ferro- 
electric materials is the high sensitivity of ferroelectric properties to the volume, 
combined with the several percent errors in predicted LDA equilibrium volumes. 
One of our main goals is to determine whether the WDA can repair this problem 
while retaining the good features of the LDA. We begin by calculating the equi- 
librium volumes of several cubic or near cubic oxides, namely CaO, SrO, BaO, 
BaTi0 3 , KNb0 3 and KTa0 3 . Some results for KNb0 3 were reported previously. 
[10] The LDA underestimates the cubic lattice parameters of these materials by 
1-2%, while at least in BaO and BaTi0 3 GGA calculations seriously overestimate 
the volume. [9] 

The calculations were done using a planewave basis set pseudopotential method, 
[31,32] with hard Troullier-Martins pseudopotentials [33] and the implementation 
of the WDA discussed in Ref. [23]. The Perdew-Wang form of G, i.e. that from 
the uniform electron gas coupling constant averaged pair correlation function was 



TABLE 1. LDA, WDA and experimental lattice parameters in A for some cubic materials. 
Values for the elements C, Si, Mo and V and are from Ref. [23]. BaTiOs and KNbC-3 are really 
rhombohedral at low temperature, but the rhombohcdral strain is small. Calculations were done 
for the cubic structure. Numbers in parentheses are the percentage deviation from experiment. 



Material 


LDA 


WDA 


Expt. 


Material 


LDA 


WDA 


Expt. 


C 


3.53 (-1.1) 


3.56 (-0.3) 


3.57 


Si 


5.36 (-1.3) 


5.40 (-0.6) 


5.43 


Mo 


3.11 (-1.3) 


3.14 (-0.3) 


3.15 


V 


2.93 (-3.0) 


2.99 (-1.0) 


3.02 


CaO 


4.71 (-2.0) 


4.81 (-0.0) 


4.81 


SrO 


5.06 (-1.8) 


5.16 (-0.0) 


5.16 


BaO 


5.46 (-1.4) 


5.56 (+0.3) 


5.54 


BaTi0 3 


3.95 (-1.2) 


4.00 (-0.1) 


4.00 


KNb0 3 


3.96 (-1.4) 


4.02 (+0.1) 


4.02 


KTaOa 


3.92 (-1.6) 


3.98 (-0.2) 


3.98 



used for the WDA. The pseudo-potentials included the semi-core states as valence 
states. In particular, 3s and 3p states of K, Ti and Ca, 4s and Ap states of Nb 
and Sr and 5s and hp states of Ba and Ta were treated in the electronic structure 
calculations, while lower lying states were pseudized. Well converged basis sets, 
including planewaves up to 121 Ry, were used, except for CaO where a higher 
cut-off of 132 Ry was used. 

One complication in the WDA is the need to use shell partitioning to avoid un- 
physical exchange-correlation interactions between valence and core electrons. It 
is used to prevent the core electrons from contributing to the exchange correlation 
hole seen by valence electrons, since it is unreasonable for tightly bound core elec- 
trons to dynamically screen valence electrons. In shell partitioning, the interaction 
between the valence electrons is treated with the WDA, but core-core and core- 
valence interactions are treated within the LDA as discussed in Refs. [23,10]. In the 
present calculations, the semi-core states on the metal ions and the O 2s states are 
treated with the LDA, lower states are pseudized with an LDA pseudopotential, 
and higher states are treated with the WDA. 

The LDA and WDA lattice parameters of these oxides are summarized in Table 
1 along with some previous results for elements. As may be seen, the WDA lattice 
parameters are quite dramatically improved relative to the LDA, including BaTi03 
and BaO, for which GGA calculations yield overcorrected volumes. 

Next we turn to phonons and the ferroelectric soft mode. The above results 
imply that the WDA can generally correct LDA errors in the equilibrium volume. 
The question then arises as to the extent that the WDA is able to reproduce the 
desirable features of the LDA with the volume corrected. LDA calculations at 
the experimental lattice parameter have demonstrated very good agreement with 
experimental data for phonons, crystal structures and even derived quantities like 
transition temperatures. Moreover, in ferroelectric/piezoelectric materials these 
quantities are governed by delicate balances between competing terms. In KNb03, 
very different results are obtained for the ferroelectric soft mode if the lattice param- 
eter is varied by 1-2 percent, and in fact the WDA does change the LDA prediction 
of the lattice parameter by an amount of this order. Thus one is led to ask whether 
the WDA will drastically change LDA results for phonons or the ferroelectric mode 



TABLE 2. LDA and WD A T-point TO phonon frequencies for cubic perovskitc structure 
KNDO3. The lattice parameter is fixed at its experimental value, which is the same as the 
WDA value. 



T-point mode 


LDA frequencies (cm r ) 


WDA frequencies (cm 1 ) 


r '^(unstable) 


195 i 


176 i 


Tl5 


170 


143 


Tl5 


473 


497 


r 25 


247 


257 



or otherwise degrade the description of the material. 

To begin addressing this issue we performed frozen phonon calculations of the 
zone center transverse modes in KNDO3. Parallel LDA and WDA planewave pseu- 
dopotential were performed as described above using a 6x6x6 special k-point sam- 
pling of the Brillouin zone. 

Calculations of atomic forces were performed for small (less than 0.05A) displace- 
ments of the atoms from the ideal cubic sites in a rhombohedral symmetry consis- 
tent with the Ti5 modes. A total of seven such force calculations were performed. 
The T 15 dynamical matrix was then least squares fit to these force calculations, 
and diagonalized to obtain the frequencies and eigenvectors. The LDA and WDA 
eigenvectors were very similar to each other and to previously reported [4,10] lin- 
earized augmented planewave (LAPW) [34] results, and are not reproduced here. 
There is only a single T 2 5 frequency, which was calculated separately, using force 
calculations with O displacements according to this mode. 

Calculated LDA and WDA frequencies of the Ti§ and r 2 5 modes in cubic per- 
ovskite structure KNbOs are given in Table 2. The LDA results are similar to 
those obtained previously by all electron LAPW, [4] pseudopotential linear re- 
sponse LAPW methods [7,8] and planewave pseudopotential calculations for the 
Ti5 modes. [35] The agreement of the LDA results with the linear response calcu- 
lations of Yu and co-workers et al. [7,8] is almost exact, with a maximum deviation 
of 4 cm -1 , while somewhat larger differences from the earlier all-electron results of 
Singh and Boyer [4] are present. Tests [3,9] have shown that most of this latter 
difference is due to the use of a 4x4x4 k-point mesh in the calculations of Ref. 
[4] leading to an underestimation of the ferroelectric instability. Fontana et al. 
[36] have performed measurements of phonon frequencies in KNbOs. These should 
not be uncritically compared to the present ground state calculations because the 
KNb0 3 undergoes three structural transitions as the temperature is raised from 
OK, and only becomes cubic at approximately 700K. Nonetheless, as discussed in 
Ref. [4], they are consistent with the LDA results for the stable modes. 

The WDA phonon frequencies are very similar to the LDA results. The main 
difference is that the splitting between the ferroelectric soft mode and the middle 
(lowest stable) r 15 mode is smaller in the WDA. However, differences between 
the LDA and WDA frequencies are comparable to the numerical accuracy of the 
calculations as estimated from the spread between the recent results of various 



groups. In any case, the present results do not provide any indication that the WDA 
degrades the favorable agreement of LDA phonon frequencies with experiment in 
these materials although further tests are clearly in order. 




FIGURE 1. Total energy as a function of ferroelectric displacement in KNDO3 within the WDA, 
LDA and GGA using both planewave pseudopotential (LDA and WDA) and LAPW approaches 
(LDA and GGA). The displacement 5 is in percent of the experimental displacement of Ref. [37]. 
Calculations are at the experimental volume, which effectively coincides with the WDA and GGA 
volumes. 



Next we turn to the ferroelectric instability. The present LDA and WDA calcu- 
lations done with a planewave pseudopotential method, and earlier LAPW calcu- 
lations both within the LDA and with GGA density functionals are in agreement 
with the experimental distortion [37] regarding the soft mode eigenvector. However, 
both LDA and GGA calculations give a soft mode amplitude that is 20% smaller 
than the existing experimental value, which was determined by powder neutron 
diffraction. Interestingly, both LDA and GGA calculations [9,10] for BaTiOs give 
a distortion that agrees very closely, (within 5%) with experiment. Fig. 1 shows 
the result of WDA calculations of this instability for KNb03 at the experimen- 
tal volume, in terms of the reported experimental distortion. As may be noted, 
there is some scatter between the curves, which represent three different exchange 
correlation potentials and two different band structure methods. However, all the 
calculations are in agreement that there is an instability of order 2 mRy/f.u. and 



that the distortion is only 80% of the reported experimental distortion. In partic- 
ular, the WDA does not significantly change the distortion relative the the LDA 
value. Taken together the results strongly suggest an experimental re-examination 
of the low temperature structure of KNbOs. 

IV LINEAR RESPONSE IN THE WDA 

As mentioned, LDA systematically overestimates the static dielectric constant, 
5oo. Not only is this itself an important characteristic of ferroelectric materials, it 
also underscores the fact that the dielectric response in general is not accurately 
reproduced, which could lead to difficulties calculating small energy differences 
associated with phonons and lattice distortions. The problem has been discussed 
in terms of the differences between the real one-electron excitation spectrum and the 
LDA band structure, primarily to underestimation of the band gap. An apparent 
paradox here is that the exact DFT theory should reproduce the static response 
functions, but not the excitation spectrum. The solution of this paradox is that 
the RPA (random phase) dielectric function, which is directly defined by the one- 
electron spectrum, is enhanced by the so-called exchange-correlation local field 
corrections. Correspondingly, a small gap and small local field corrections result 
in the large gap and large corrections. These corrections, in turn, 

are defined by the second variation of the exchange-correlation energy, K xc (r, r') = 
S 2 E xc /Sn(r)Sn(r').The LDA provides correct K xc (r,r') only in metals and only in 
the long range limit. 

One of the first attempts to correct this was the formulation of ADA, where 
n(r') in Eq. (1) is substituted by n(r), so that E A ® A = J n(r)e xc [n(r)]dr. Then 
n(r) is defined as n(r) = / w[\r — r'|, h(r)]n(r')dr', and the universal function w is 
chosen so that 5 2 E£P A /5n(r)5n(r') gives the correct K xc (r, r') for the uniform gas. 
Contrary to the WDA, the ADA is not self-interaction free in one electron systems, 
and thus was never as popular as WDA. 

From the beginning there was substantial interest in the behavior of WDA in 
the delocalized limit [27]. Williams and von Barth [38] suggested that the WDA 
should give substantial improvement over the LDA in this limit, but till now no 
systematic study has been reported. If this conjecture is true, the WDA has a great 
advantage over any other known approximation to the DFT in the sense that it 
would accurately reproduce two key physical limits. At least some improvement 
over LDA is to be expected: in the short range K^ DA (r, r') remains finite, and 
correspondingly decay with q — > oo in reciprocal space. Smaller K xc will result 
in weaker local field corrections and in an improvement in e^. It is not clear a 
priori, though, how much K x v c DA { y q) is improved over the LDA over the whole q 
range and, if the improvement is only modest, whether or not an approximation 
based on the WDA exists that does provide proper behavior. Here we derive an 
expression for K xc in the WDA, calculate K xc for popular flavors of WDA, and 
discuss construction of a WDA method with improved K xc . 



We start by deriving a closed expression for K xc in the WDA for arbitrary G. First 
some notation: denote the product (e 2 /r)G(r) as W(r), use atomic units (e = 1, 
h — 1), and use primes for the derivative with respect to the density argument, 
e.g. G' = dG/dn. We also introduce two functions, reflecting implicit dependence 
of the weighted density n on variations of the real density: 

d(r'-r) = 8n{v')/8n{v) (3) 

/(r'-r r'-r") = ^ = (4) 
M ' ; 5n(r)5n(r») 8n{v") ' K ' 

Using the WDA expression for the exchange-correlation energy, 

E xc = (1/2) J n(r)n(r')W[\r-r'\,n(r)]drdr', (5) 

we express K xc in terms of functions d and /, and find these functions using the 
normalization condition. 

J dr'n(r')G[\r — r'|, n(r)] — — 1. (6) 

We proceed then in reciprocal space, which corresponds to using density perturba- 
tion of the form Sn(r) =n q e iqr . Let W q , G q , d q and f Ptq be the Fourier transforms 
of the corresponding functions. Then 

dg = -G q /ng' . (7) 

Since at q — > the LDA should be restored, 

dr'W[\r-r'\,n]=2e xc /n. (8) 



From this it immediately follows that 

G = -I /n, W = 2e xc /n. (9) 

Thus d q = —nG q . Next variation of Eq.(6) gives us f Piq . In fact, we need only 
diagonal elements, f q - q , for which we find f q - q = 2nG q {nG' q + G q ). The second 
variation of Eq. (5) in terms of d and / is 

17? 

K xc (q) = W q + nd q W' q + nd q W' Q + y (« + n 2 /,,-X)- 

resulting in 

K xc (q) =W q - n 2 G q {W' q + + n 2 (n 2 G 2 ^)'/2. (10) 

The original formulation of the WDA used the homogeneous electron gas function 
for G. Since then, three forms of G have been used in calculations, all of which result 



in improvement over LDA (in the admittedly limited number of tests performed to 
date). These are: the function G derived for the uniform gas by Perdew and Wang 

[28], the Gunnarsson- Jones function G GJ (r) = Ci(ra){l — exp[— (^^y) ]}, and 

the Gritsenko et al. [39] function G GRBA (r) = d(n) exp[- (^y)^}, k = 1.5 (note 
that the uniform gas function [28] is approximately given by the same expression 
with k = 2). We tested these functions for the densities r s = 1, 2, and 5 and 
obtained modest agreement with the Monte Carlo results [40] (Cf. the left panel of 
Fig. 2, where we plot the calculated exchange-correlation local field factor I xc (q) = 

2 

j^K xc (q), and compare it with Monte Carlo data [40]). By construction, K xc (0) is 
correct (and is in fact the LDA value). At q « 1.5 — 1.8k F K xc falls below its LDA 
value and continues to decrease at large g's. However, a closer look reveals some 
disagreements: first, I^ DA (q) is considerably larger than the Monte-Carlo data for 
the wave vectors between « 0.5k F and 1.5k F . Second, I xc (q) in WDA tends to a 
constant value, while in Monte Carlo calculations it is K xc (q) itself that has a finite 
limit at q — > oo, and I xc (q) — > const ■ q 2 at q — > oo. 

Can one correct these two deficiencies without compromising the correct one- 
electron limit of WDA? As discussed above, there is no particular reason to use 
the homogeneous electron gas pair correlation function for G (nor, as discussed 
above, the exact pair correlation function for the inhomogeneous system, even if it 
had been known). Since using Go in WDA does not guarantee any improvement 
in describing properties of the homogeneous gas itself, one may use the freedom in 
G(r) to adjust the WDA so that the calculated local field factor (and thus linear 
response function) is as accurate as possible. Inversion of eq.(10) yields G(q) for 
a given K xc (q). It does not guarantee, however, that the result will be physical. 
So, as a first step, let us analyze Eq. (10). Let us first mention that the real 
K xc {q) changes its behavior from the long range limit to the short range limit near 
q = 2k F , which plays the role of the inverse length scale. What is the characteristic 
length scale in WDA? To find that, we write G q = —<p(p/Q)/n, with the condition 
y?(0) = 1, where Q is some constant (both the Gunnarsson- Jones and the Gritsenko 
et al. functions are of this form). Then 

W = 1 r G q dq = r <p( X )d X = *=*5>. (11) 

7r Jo ' nn Jo n 

If we now define Q(n) = —7ie xc (n), then the second condition on ip(x) becomes 
/ °° ip(x)dx = 1. These two conditions reduce our freedom to adjust G q : since 
the characteristic size of (p(x) is of order of 1, the wave vector dependence of G q 
is defined by the ratio q/Q = —q/iie xc . Apparantely, varying the shape of the 
function (p will not significatntly change the lentgh scale of the resulting K xc (q). 
Explicitly density-dependent functions <p(x,n) may be needed to shift the hump 
from its position of q ~ 1.5k F to q ~ 2k F . It is still an open question whether or 
not a physically sound function can be found with this property. 

However, even if the u 2k F " problem is fixed, another, probably even more im- 
portant problem remains: the short wave length behavior of K xc . It is easy to see 



that if G(q) — > at q — > oo, then H^, — >■ const/p 2 at p — > oo, and so does, according 
to Eq.(lO), K xc . On the other hand, as mentioned above, the correct K xc (q) goes 
to a constant at q — > oo as g 2 , although the constant is smaller than K^ A . This 
result was predicted by Holas [41] and is physically important: it comes from the 
exchange-correlation contribution to kinetic energy (which is essentially local and 
decays slower with q than the interaction part of E xc ). The present WDA misses 




FIGURE 2. Left panel: exchange-correlation local field factor in the WDA of Ref. [30] (Gun- 
narsson- Jones), Ref. [39] (Gritscnko ct al.), and derived from the homogeneous electron gas pair 
correlation function (Pcrdcw-Wang) , as compared with the Monte Carlo results (Monte Carlo) 
and the interpolating formula thereof (MC interpolation), as given in Ref. [40]. Densities, from 
top to bottom, correspond to r s = 1, 2, 5. Right panel: the same for the modified WDA of Eq.16. 
Also the analytical formula of Farid et al. Ref. [42] is shown. 



the corresponding physics. Fortunately, this is easy to correct. Farid et al. [42] 
tabulated the coefficient 7 that defines the asymptotic behavior of K xc (q — > 00) 
as K xc (q — > 00) = ~^l( n )§r, where 7(71) is a universal function, parametrized in 
Ref. [42]. Let us now modify the function G(r) 

G{r) = Gi(r) + G 2 (r) = A5(r)/4irr + G 2 (r), (12) 

Since / G\{r)r 2 dr = 0, the normalization condition for G 2 is the same as for G 
itself. Since 4n J Gi{r)rdr = A, the LDA limit condition for G 2 becomes 

4tt J G 2 (r)rdr = 2e xc (n)/n, 

e xc (n) = e xc (n) - An/2. (13) 

Thus 

A = -47T7(n)//4, 
txc{n) = e xc (n) + 2nn-f(n)/k 2 F 

Now G p = G 2p , W P = A + W 2p , and W p = A' + W' 2p , 

K xc {q) = A + W 2q - n 2 G 2q (A' + W' 2q + W^) + n 2 (n 2 G^ 2 ' )' 

= A - n\G 2p A' + (15) 

where K xc is calculated from e xc in exactly the same way as K xc is calculated from 
e xc . The corresponding functional for the exchange-correlation energy is 

E^° A = I J ^^G[\v - r'\,n(r)]drdr' + / n(r)^[n(r)]dr. (16) 

Here G(r) is normalized to e :rc (n) instead of e xc (n). In practice, the first term gives 
rise to the standard expression for the WDA potential [22,21], and the second yields 
two additional terms, one from the variation of n(r), and the other arising from 
5n(r')/Sn(r). Since we do not require that G(r) = fo[g(r; A, n) — l]d\., where 
g corresponds to the uniform gas, but rather consider it to be a flexible function 
satisfying two normalization conditions, further improvement of the method should 
be possible along the line described in the previous paragraph, namely the freedom 
in choosing G(r) can be used to yield K xc according to Eq. (15) close to the linear 
response of the homogeneous electron gas, including correct behavior near q = 2kp- 
In the right panel of Fig. 2, we show I xc calculated according to Eq. 16 with the 
different functional form of G(r). Clearly, the results are much better than either 
the LDA or "conventional" WDA. 

In short, we have calculated the exchange-correlation local field function K xc in 
the WDA, and found that besides the expected improvement over the LDA it has 
two major deficiencies: (1) it does not have correct asymptotic behavior at q — > 00, 



(14) 



and (2) the characteristic feature at q rs 2k F is displaced towards smaller g's. 
The former can be easily corrected by adding a delta-function component to G(r), 
which results in Eq. (16). The latter is harder to fix, but there are still unused 
degrees of freedom in the formalism which may be used to tune the behavior near 
2kp- Intuitively (cf. Ref. [38]), a method which retains exact one-electron limit 
of WDA, and at the same time is accurate in the opposite limit of the nearly 
uniform electron gas, seems promising for practical applications. However, tests on 
real materials will be needed to determine whether or not this modification of the 
WDA is advantageous in practice. 

V CONCLUSIONS 

WDA calculations of the equilibrium volume of several oxides show that with 
the electron gas form of G and shell partitioning, the WDA yields much improved 
volumes over the LDA and GGA. No case was found where the WDA degrades 
the LDA results. Phonon frequencies and the ferroelectric distortion in KNbO"3 
are in good accord with LDA predictions provided that the LDA calculations are 
performed at the experimental volume, which is effectively the same as the WDA 
volume. The direction of the ferroelectric soft mode is in good agreement with 
experiment, but all exchange correlation functionals tested yield a distortion that 
is 20% smaller than that reported in the powder neutron experiment of Hewat. [37] 
The implication is that this displacement should be re-examined from an experi- 
mentally. The above results combined with the fact that the WDA linear response 
of the electron gas is in substantially better agreement with Monte Carlo data than 
the LDA, and that, if desired, the approach has the flexibility to further improve 
this property, is suggestive that the WDA may be a generally more reliable method 
than the LDA. 

We are thankful for helpful discussions with L.L. Boyer, H. Krakauer, O. Gun- 
narsson and R. Resta. This work is supported by the Office of Naval Research. 
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